# Paper: Is it still the economy? Economic voting in polarized politics
# Individual-level simulation
# Replication: Figure S7 (Supplementary Material)
remove(list = ls())

library(conflicted)
library(tidyverse)
library(dplyr)
conflict_prefer("filter","dplyr")
library(ggplot2)
library(here) # When we have an R-project, it reads files without setting an directory
conflict_prefer("here", "here")
library(nnet)
library(ggeffects)
library(MASS)
conflict_prefer("select","dplyr")

load(here("Final_Paper", "R&R", "IndLevel.RData"))

repinc_data <- model_data %>% filter(rep_inc == 1)
table(repinc_data$year)
deminc_data <- model_data %>% filter(rep_inc == 0)
table(deminc_data$year)

table(model_data$econ_eval)

table(repinc_data$state_fips)

repinc_data2 <- repinc_data %>% 
  filter(!is.na(econ_eval)) %>% 
  filter(!is.na(percep_polar)) %>% 
  filter(!is.na(ptd_id3)) %>% 
  filter(!is.na(female)) %>% 
  filter(!is.na(black)) %>% 
  filter(!is.na(hispanic)) %>% 
  filter(!is.na(asian)) %>% 
  filter(!is.na(other_race)) %>% 
  filter(!is.na(log_age)) %>% 
  filter(!is.na(college)) %>% 
  filter(!is.na(income)) %>% 
  filter(!is.na(state_fips)) %>% 
  mutate(elec_04 = ifelse(year == 2004, 1, 0)) %>% 
  mutate(elec_08 = ifelse(year == 2008, 1, 0)) %>% 
  mutate(elec_20 = ifelse(year == 2020, 1, 0))

table(repinc_data$year)
table(repinc_data2$year)
summary(ml_repinc <- MASS::polr(econ_eval ~ percep_polar*ptd_id3 + 
                                  female + black + hispanic + asian + other_race + 
                                  log_age + college + income + 
                                  elec_04 + elec_08 + elec_20, 
                                data = repinc_data2, Hess=TRUE))
ml_repinc


deminc_data2 <- deminc_data %>% 
  filter(!is.na(econ_eval)) %>% 
  filter(!is.na(percep_polar)) %>% 
  filter(!is.na(ptd_id3)) %>% 
  filter(!is.na(female)) %>% 
  filter(!is.na(black)) %>% 
  filter(!is.na(hispanic)) %>% 
  filter(!is.na(asian)) %>% 
  filter(!is.na(other_race)) %>% 
  filter(!is.na(log_age)) %>% 
  filter(!is.na(college)) %>% 
  filter(!is.na(income)) %>% 
  filter(!is.na(state_fips)) %>% 
  mutate(elec_00 = ifelse(year == 2000, 1, 0)) %>% 
  mutate(elec_12 = ifelse(year == 2012, 1, 0)) %>% 
  mutate(elec_16 = ifelse(year == 2016, 1, 0))

glimpse(deminc_data2)
table(deminc_data2$year)
table(deminc_data$year)
table(deminc_data$state_fips)


summary(ml_deminc <- MASS::polr(econ_eval ~ percep_polar*ptd_id3 + 
                                  female + black + hispanic + asian + other_race + 
                                  log_age + college + income + 
                                  elec_00 + elec_12 + elec_16, 
                                data = deminc_data2, Hess=TRUE))

ml_deminc


ml_repr <- ml_repinc %>% 
  ggeffect(terms = c("ptd_id3", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_repr)
ml_repr <- data.frame(ml_repr)

ml_repr <- ml_repr %>% 
  mutate(incumbent = "Rep")


ml_demd <- ml_deminc %>% 
  ggeffect(terms = c("ptd_id3", "percep_polar[0.1, 0.5, 1]"), 
           ci.lvl = 0.95)
glimpse(ml_demd)
ml_demd <- data.frame(ml_demd)

ml_demd <- ml_demd %>% 
  mutate(incumbent = "Dem")

ml_inc <- rbind(ml_repr, ml_demd)
glimpse(ml_inc)


ml_inc <- ml_inc %>% 
  rename(part_id = x)

ml_inc <- ml_inc %>% 
  mutate(econ_eva = ifelse(response.level == "X1", "Worse", NA),
         econ_eva = ifelse(response.level == "X2", "Same", econ_eva),
         econ_eva = ifelse(response.level == "X3", "Better", econ_eva))

g_rep <- ml_inc %>% 
  filter(incumbent == "Rep") %>% 
  mutate(part_id = factor(part_id,
                          levels = c("Dem", 
                                     "Ind",
                                     "Rep"),
                          labels = c("Dem", 
                                     "Ind",
                                     "Rep"))) %>% 
  mutate(econ_eva = factor(econ_eva,
                           levels = c("Worse", 
                                      "Same",
                                      "Better"),
                           labels = c("Worse", 
                                      "Same",
                                      "Better"))) %>% 
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(econ_eva == "Worse") %>% 
  #mutate(x = as.factor(x)) %>%
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(colour = part_id, shape = group, fill = part_id)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = part_id, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .65, position = position_dodge(width = .5)) + 
  theme_light() + 
  labs(x = NULL, y = "Predicted Probability", 
       title = "(a) Republican Presidency") +
  scale_color_manual(name = "Partisanship:", 
                     values = c("blue", "gray37", "firebrick3"),
                     labels = c("Dem", "Ind", "Rep")) +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", "0.5", "1"),
                     values = c(22, 21, 24)) +
  scale_fill_manual(name = "Partisanship:",
                    values = c("blue", "gray37", "firebrick3"), 
                    labels = c("Dem", "Ind", "Rep")) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x=element_blank(),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 11),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))
g_rep

g_dem <- ml_inc %>% 
  filter(incumbent == "Dem") %>% 
  mutate(part_id = factor(part_id,
                          levels = c("Dem", 
                                     "Ind",
                                     "Rep"),
                          labels = c("Dem", 
                                     "Ind",
                                     "Rep"))) %>% 
  mutate(econ_eva = factor(econ_eva,
                           levels = c("Worse", 
                                      "Same",
                                      "Better"),
                           labels = c("Worse", 
                                      "Same",
                                      "Better"))) %>% 
  mutate(group = as.factor(group)) %>% 
  #mutate(response.level = as.factor(response.level)) %>%
  filter(econ_eva == "Worse") %>% 
  #mutate(x = as.factor(x)) %>%
  #filter(response.level == "Inc") %>% 
  #filter(group == 0.33 | group == 1) %>% 
  ggplot(aes(colour = part_id, shape = group, fill = part_id)) + 
  #geom_hline(yintercept = 0, colour = "black", lty = 2, size = .75) +
  geom_pointrange(aes(x = part_id, y = predicted, 
                      ymin = conf.low, ymax = conf.high),
                  lwd = .65, position = position_dodge(width = .5)) + 
  theme_light() + 
  labs(x = NULL, y = NULL, 
       title = "(b) Democratic Presidency") +
  scale_color_manual(name = "Partisanship:", 
                     values = c("blue", "gray37", "firebrick3"),
                     labels = c("Dem", "Ind", "Rep")) +
  scale_shape_manual(name = "Ideological gap:",
                     labels = c("0.1", "0.5", "1"),
                     values = c(22, 21, 24)) +
  scale_fill_manual(name = "Partisanship:",
                    values = c("blue", "gray37", "firebrick3"), 
                    labels = c("Dem", "Ind", "Rep")) +
  scale_y_continuous(breaks = c(1, 0.8, 0.6, 0.4, 0.2, 0),
                     limits = c(0,1)) +
  theme(axis.text.x=element_blank(),
        axis.text.y = element_text(size = 10),
        axis.title.y = element_text(size = 12),
        legend.text = element_text(size = 11),
        legend.title = element_text(size = 11),
        legend.position = "bottom",
        legend.direction = "horizontal",
        plot.title = element_text(size = 14))
g_dem


multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

multiplot(g_rep, g_dem, cols = 2)

